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Stationary points (SPs) of the potential energy landscapes can be classified by their Morse index, i.e., the 
number of negative eigenvalues of the Hessian evaluated at the SPs. In understanding chemical clusters 
through their potential energy landscapes, only SPs of a particular Morse index are needed. We propose a 
modification of the “fixed-point homotopy” method which can be used to directly target stationary points of a 
specified Morse index. We demonstrate the effectiveness of our approach by applying it to the Lennard-Jones 
clusters. 


I. INTRODUCTION 

The stationary points (SPs) of the potential energy 
function, U(x), with x = (xi,...,x n ), are the bases 
of the potential energy landscape methods that have 
recently gained significant attention from the chemical 
physics community 1 . Finding SPs of U(x) amounts to 
solve the corresponding system of nonlinear equations 
= 0, i = 1 The stability and the local 

geometry of the SPs are determined by the eigenvalues 
of the Hessian matrix of V (x), %(U)(x)ij = d 2 V/dxidxj, 
evaluated at the SPs. While the eigenvalues themselves 
depend on the choice of coordinate systems, by the 
Sylvester’s Law of Inertia, their signs are independent 
from coordinate systems and hence are intrinsic to the 
local geometry of the SPs. The number of negative eigen¬ 
values of %(U)(x) is known as the Morse index or simply 
index. An SP x is said to be degenerate if FI(U)(x) is 
singular and nondegenerate otherwise. 

A plenty of numerical methods to solve such systems 
of nonlinear equations® ® as w ell as methods to certify 
numerical approximation^ 9 20 ^are available in the energy 
landscape areas and beyond (see Ref. [2T)for a recent short 
survey). Among these, homotopy continuation methods 
have long been used to solve multivariate nonlinear sys¬ 
tems of equations (see Ref.[22lfor surveys in this subject). 
Among a rich body of works applying homotopy meth¬ 
ods to scientific problems (e.g. Ref. 123$, recently, the nu¬ 
merical polynomial homotopy continuation methocf®^ 31; 
has gained special attention for its ability to find all the 
complex and real SPs of the potentials with polynomial- 
like nonlinearities. However, when one is interested in 
finding only real solutions, this method may turn out to 
be computationally expensive. The authors have applied 
another homotopy based method, the Newton homotopy, 
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which can be rest ricted to find only real solutions, with 
noteworthy result J® 2 33 . However, this method, same as 
almost all homotopy based methods, does not distin¬ 
guish among SPs of different indices whereas in chemical 
physics applications, finding stable structures and other 
thermodynamic properties of the systems utilizes only 
on SPs of indices zero (called minima) and one (called 
transition states^. 

Among the great variety of homotopy methods, the 
fixed-point homotopy method®®, being one of the first 
general homotopy methods invented, has long been sus¬ 
pected to be able to target index 0.® Modifications to the 
fixed-point homotopy method have been proposed to tar¬ 
get a wider range of SPs, especially index 1 SPs. Most 
notable among them was the “singular fixed-point ho¬ 
motopy” method®] which uses index 0 SPs as bifurcated 
starting points and is likely to be able to reach nearby 
SPs of index 1. 

Several other methods have been developed which at¬ 
tempt to find SPs of specific index. Some of the more di¬ 
rect methods include eigenvector following methods 1 -®®, 
single- and dou ble-ende d searches and gradient square 
based method jH l 10 l 12 l 39 ], e tc. However, most of these 
methods are only locally convergent, i.e., one needs to 
start from an initial guess close to the actual SP. More¬ 
over, some of these methods also require at least one SP 
to start with. 

In the present contribution, we propose a novel mod¬ 
ification to the fixed-point homotopy which can directly 
target SPs of a given index. Unlike the “singular fixed- 
point homotopy” method, our proposed method starts 
with random starting points and does not require an ex¬ 
isting collection of index zero SPs as starting points (cf. 
Ref. 1231) . and it is hence more likely to provide a global 
sample of index 1 SPs. This feature is of particular im¬ 
portance when no index zero SP is known. Compared to 
other methods based on gradient flow or quasi-Newton 
techniques, the proposed method is also likely t o hav e a 
strong advantage in dealing with degenerate SP d 321331 . 
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II. THE DEFORMATION OF AN ENERGY LANDSCAPE 

In the following, V (x) the potential enrgy function to be 
explored will be called the target potential energy func¬ 
tion. For an integer m with 0 < m < n, called the target 
index , our goal is to locate SPs x = (aq,... ,x n ) £ R” 
with index exactly to. That is, we would like to find 
points x £ R” such that 

1. W(x) = 0; and 

2. H(V) at x has m negative eigenvalues. 


above, we construct a family of closely related potential 
V t parametrized by t so that V t = V yet the SPs of the 
starting potential Vo can be located and studied easily. 
Here the starting potential is chosen to be the “standard 
?n-saddle” located at a point a = (ai,... ,a n ) which is 
represented by the potential 

-(ix-ap 2 -- (x m -a m ) 2 + (x rn + 1 -a m+1 ) 2 + ■ • - + (x„-a„) 2 . (1) 

Clearly, x = a is a nondegenerate SP of the above func¬ 
tion with index exactly to. For brevity, we define the 
quadratic form as 


Criterion 1 amounts to solving a system of non¬ 
linear equations. In the present article, we focus 
on a h omotopy con tinuation approach for solving this 

problenJ^MIEMlllol 

in which the target potential func¬ 
tion V is embedded into a family of closely related po¬ 
tential functions V t parametrized by a parameter t so that 
at t = 1, Vi = V. Geometrically speaking, this can be 
understood as the continuous deformation of the target 
energy landscape among a family of related landscapes. 
The deformation is constructed so that the SPs of one 
member Vo at t = 0 can be studied easily. The move¬ 
ment of the SPs of Vo under this deformation as t varies 
from 0 to 1 can then be tracked. In particular, robust 
numerical algorithms can be used to trace the trajectory 
of the SPs. If the trajectory can be extended to t = 1, 
an SP of the target potential is then produced as V\ = V. 

In certain cases, the index can be preserved along the 
trajectory formed by the aforementioned deformation. 
By the continuity of eigenvalues of 'H(V t ) along a tra¬ 
jectory it is easy to see that the index must remain the 
same unless a degenerate SP of V t is encountered where 
at least one eigenvalue of H(V t ) vanishes^ However, this 
criterion is quite limited, since from the point of view of 
catastrophe theon ^2143] ag one traces the SPs 0 f a f am ily 
of potential functions, an encounter with a degenerate 
SP is not only possible, but sometimes necessary. 

In the following we propose a new homotopy construc¬ 
tion that can potentially target SPs of a given index un¬ 
der much relaxed conditions. 


III. INDEX-RESOLVED FIXED-POINT HOMOTOPY 

The fixed-point homotopj^USSl j s one 0 f the first ho¬ 
motopy methods devised to solve general nonlinear sys¬ 
tems. Here we propose an 11 index-resolved fixed-point ho¬ 
motopy” method for finding SPs of a given index. As 
we shall demonstrate in §IV| with experiments on the 
Lennard-Jones cluster, it is capable of quickly obtaining 
a large number of SPs of a given index. 

In keeping with the general framework for finding SPs 
of a smooth potential using homotopy methods described 


jW{ u) 



u 1 



( 2 ) 


where u = (ui,... ,u n ), and I m and I n -m are identity 
matrices of dimension to x to and (n — to) x (n — to), 
respectively. The standard ?n-saddle 0 at x = a is then 
represented by the quadratic function j( m )(x — a). 

Representing the deformation between this m-saddle 
and the target potential V is the one-parameter family 


V* (m) ( x ) := (1 - t) J (m) (x - a)+fV(x). (3) 


This family contains the target potential V (x) as a mem¬ 
ber, since at t = 1, V^™"* = V. The starting potential 
Vj^ m ^(x) = J< m )(x — a), at t = 0, has a unique SP of 

index to. As t varies, represents a smooth deforma¬ 
tion between the two potentials. 

We now consider the effect of the deformation <[3j) on 
the SP x = a, that is, how the SP x = a of Vq" 1 ' 1 relates 

to SPs of nearby members of Vf"' 1 as t moves away from 
t = 0. The local theory of smooth real-valued functions 
yields that under a small perturbation, nondegenerate 
SPs remain nondegenerate, i.e., by varying t sufficiently 
close to 0, the nondegenerate SPs of V^ m> traces out a 
small piece of smooth curve. The homotopy continuation 
approach hinges on the “continuation” of this small piece 
of curve into a globally defined curve that could poten¬ 
tially connect the SP x = a of V^ m> to an SP of the target 
potential V^ m ' > = V which we aim to find. A standard 
application of the Generalized Sard’s Theoren guar¬ 
antees that this is possible for “generic” choices of a. 
Putting aside the technical statement of the theorem, it 
essentially means that if a, the (unique) SP of the start¬ 
ing potential, is chosen at random, then with probabil¬ 
ity on the set of all SPs of for t £ [0,1) form 
of smooth curves in the (n + 1)-dimensional x-f space. 
Among them, we focus on the curve that passes through 
x = a and t = 0, which will be called the trajectory of a. 


This is proved in Ref. 1371 and used without explicit statement in 
several related works e.g. Ref. 1341 and im 


Equivalently, this is stating that the probability of choosing the 
“bad” start points is zero. 
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Since V} m ' > = V, if this trajectory extends to t = 1 then 
an SP of the target potential is obtained. 

From the algebraic standpoint, if we define the index- 
resolved, fixed-point homotopy with index m as 


FW(x,i):=Vi7*(x), (4) 


then the crux of the method is to numerically trace 
along the smooth trajectory which is a curve defined by 
ij( m )(x, t) = 0 and containing the point (a, 0) in the 
hope that a point (x^, 1) can be reached]^] 

While there is a rich selection of numerical meth¬ 
ods for tracing such curves, within the community of 
homotopy methods, a special class of methods based 
on a predictor-corrector scheme together with an arc- 
length parametrization of the curve has emerged as 
the meth od of cho ice for its superior stability and 
efficiency 2 ^^ 42 * 40 ^ 44 ^ ^. 


Our goal, however, is to locate SPs of index m. Since 
the starting point x = a is an index-m SP of the starting 
potential Vq , the hope is that this index is inherited by 
the SP of the target potential the trajectory obtains (if 
the trajectory reaches t = 1). Clearly, if the trajectory 
encounters no degenerate SP, the continuity of eigenval¬ 
ues of (which must remain nonzero without en¬ 

countering a degenerate SP) ensures preservation of index 
along the entire trajectory, and consequently the SP of 
the target potential it reaches has index rn. However, our 
numerical experiments with the Lennard-Jones potential 
((IV) suggest that the index can be preserved even when 
the assumption in this theorem fails to hold, that is, when 
the trajectory encounters some turning points. Some ob¬ 
servations are presented in jfV| 


A. Summary of the Proposed Numerical Algorithm 


i.e., if it passes through a point (x, t) = (xO, 1), then an 
SP of the target potential function V is producccj/] 


B. An example 

As an example, consider the potential function 

V(x,y) = -x 3 +y 3 + 2xy -6, (5) 

whose contour plot is shown in Figure [I] It has two SPs: 
a local minimum (index 0) at (—2/3, 2/3) and a saddle 
point (index 1) at (0,0). 


0.8 
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Figure 1. Contour plot of potential function § which has 
a local minimum (index 0) at (—2/3, 2/3) and a saddle point 
(index 1) at (0,0). Starting from the exact same point the 
two index-resolved fixed-point homotopies with different tar¬ 
get indices Vt and Vt define two different trajectories marked 
by the red/light and the blue/dark curves respectively. 



For the target potential function V (x) and a given tar¬ 
get index m < n, one first randomly pick a starting point 
a £ with which the family 

V t im \ x ) : = (!“*) Jm(x-a) + fV(x) 

is constructed. We are interested in the SPs of v} 1 " 1 for 
all f as a set in the (n + l)-dimensional space. This set is 
defined by the nonlinear system ff( m )(x,f) := = 

0. The starting point (a, 0) is clearly in this set. More¬ 
over, for generic choices of a, this system defines smooth 
trajectory with the starting point x = a and t = 0. 

One can then apply efficient numerical methods based 
on the predictor-corrector scheme to trace along the 
smooth trajectory. If the trajectory extends to t = 1, 


! Though not limited to gradient systems, the original fixed-point 
homotopjHHH] can be considered as a special case of this con¬ 
struction with the target index 0, i.e., the homotopy 


(Index 0) To locate the index 0 SP we construct 

Vt 0) : = —yy—[(£ — ai) 2 + (y — CI 2 ) 2 ] + tV(x,y) (6) 

according to © with target index m = 0 where 
(ai,a 2 ) = (—0.494332,0.804689) is a randomly chosen 


4 To control the “numerical condition” of the curve tracing algo¬ 
rithm, this formulation of the deformation V .can be general¬ 
ized to 

y/ m) (x) := (1 - t) jW(P(t) • (x - a)) +tvm ) ■ x) 

where for each t value P(t ) and Q(t) are nonsingular n x n ma¬ 
trices which can be chosen to tame the numerical condition of 
the corresponding system of equations VV^ m ^(x) = 0. For ex¬ 
ample, by choosing P{t) and Q(t) to be diagonal matrices with 
positive entries, one can selectively scale the individual variables 
xi,..., x n to ease the potential problems of imbalance in mag¬ 
nitudes. 
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starting point. Clearly, (x, y ) = (ai, a 2 ) is an SP of index 
0 of the starting potential Vq 0 ^ = (x — ai) 2 + (y — CZ 2 ) 2 . 

The projection of the smooth trajectory, defined by 
H^°\x,y) := VVt°\x,y) = (0,0), onto the x-y plane 
(where t is removed) is shown in Figure[l]as the red/light 
curve. It converges to the index 0 SP at (—2/3,2/3). 

(Index 1) To locate the index 1 SP, we construct 

^t (1) : = ~^~[—(^ — a i ) 2 + (y ~ a 2) 2 ] + tV(x,y) ( 7 ) 

with the only difference from ([6]) being the sign of 
(x — ai) 2 . With the same starting point ( 01 , 02 ) = 
(—0.494332,0.804689), H^(x,y) := VV t (1) (x, y) = (0,0) 
defines a different trajectory whose projection onto the 
x-y-plane is shown in Figure [l] as the blue/dark curve. 
This curve converges to the saddle point (0,0) of V(x, y) 
instead. It is worth pointing out that even though the 
starting point (01,02) is very close to the index 0 SP, as 
shown in Figure [T] the trajectory still converges to the 
index 1 SP which is much further way. 

In both cases, the target indices in in © and ([7]) de¬ 
termine which SP the corresponding curve converge to. 


IV. THE LENNARD-JONES CLUSTER 


We app ly our method to find SPs of the Lennard-Jones 
clustered of N atoms whose potential function is given as 


N N 
VN=*e± 

i =1 j=i +1 


( 8 ) 


where N is the number of atoms, 2 1 / 6 <t is the equilibrium 
pair separation, e is the pair well depth, and 


nj = \J(Xi - Xj ) 2 + (y t - y 3 ) 2 + (zi - Zj ) 2 

is the Euclidean distance between an i-th and j-tii atoms. 
We take e = <7 = 1 for simplicity. In addition to the fact 
that the Lennard-Jones cluster potential serves as a good 
approximate for the atomic interactions, the potential 
exhibits complicated landscape structure, i.e., the num¬ 
ber of minima exponentially grows when increasing N 
and has a multi-funnel structure due to the simultaneous 
presence of competing growth sequences 1 . This model 
has been under extensive searches for minima and higher 
index SP d 13 ! 48 ! which provided us bases for comparison. 

Defined in terms of the pairwise distances, Vjv is clearly 
invariant under rotation and translation. Consequently 
in the x-y-z coordinates, all SPs would have certain de¬ 
grees of freedom. After a translation of the cluster, we 
can fix the first atom at the origin, i.e., x\ = yi = Z\ = 0. 
For the ease of computation, we shall restrict our atten¬ 
tion to the cases where the N atoms are not collinear. 
We can therefore fix t /2 = Z 2 = 23 = 0 and also require 


2/3^00 After the restriction, we have a total of 3 N — 6 
variables in Vn ■ Due to the permutation symmetry, Vjy 
may exhibit the same value at different SP s. Su ch SPs 
are known as permutation-inversion isomers ! 49 ! 50 ! 


We applied the index-resolved fixed-point homotopy 
to the problem of finding SPs of the Lennard-Jones po¬ 
tential Vjv with certain indices. Since SPs of indices 0 
(local minima), 1 (transition states), and 2 are most fre¬ 
quently used in computational chemistry, we restrict our 
attention to honrotopies H (as defined in 0 ) with 
m = 0 , 1 , 2 , although general constructions with any 
to < 3 N — 6 are possible!—!—! Table 0 shows the num¬ 
ber of SPs of Vjv for a range of N values obtained by 
ij(°), H^\ and H^ respectively. Remarkably, the in¬ 
dices of all SPs found match exactly the target index 
m used in the construction of H^ m \ suggesting that the 
index-resolved fixed-point homotopy proposed here has a 
great potential in targeting SPs of a specific index. 


As a homotopy-based method, the index-resolved 
fixed-point homotopy has a strong advantage in deal¬ 
ing with degenerate SPsF^l The runs that produced SPs 
listed in Tabic Q] have also resulted in a number of SPs 
that appear to be degenerate (and hence not counted). 
Table [TT] lists a selection of these numerically degenerate 
SPs where Hessian matrices ’H(Vjv) have eigenvalues of 
magnitude less than 10~ 10 . Though the physical mean¬ 
ing of the degenerate SPs in the Lennard-Jones clusters 
appear to be understudied, their very presence highlights 
the richness of the energy landscape and merits further 
analysis. 


5 If the N atoms are not collinear, after relabeling the N atoms 
and rotations, we can always find a representative that satisfy 
these conditions. For the pathological collinear configurations 
(where all N atoms lie on a line) the restriction t /2 = -2 = ^3 = b 
does not remove all the degrees of freedom. Consequently such 
configurations are best handled by alternative formulations. 
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Target 

index 

Obtained 

index 

N 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

m = 0 

0 

1 (22) 

2 (59) 

4 (101) 

8 (134) 

15 (101) 

28 (132) 

45 (46) 

58 (193) 

68 (110) 

60 (92) 

70 (78) 

m = 1 

1 

2 (53) 

3 (92) 

8 (105) 

13 (148) 

22 (100) 

57 (391) 

65 (243) 

80 (385) 

79 (313) 

70 (96) 

64 (242) 

m = 2 

2 

4 (44) 

6 (168) 

8 (142) 

14 (170) 

14 (96) 

15 (106) 

56 (124) 

52 (143) 

57 (119) 

55 (115) 

52 (97) 


Table I. Numbers of geometric configurations of N atoms (together with numbers of SPs representing them) with different 
indices of the Lennard-Jones potential Vv (j8j for a range of N values obtained by index-resolved fixed-point homotopies with 
target indices m = 0,1, 2 respectively. E.g., the entry “2 (59)" for N = 6, m = 0, and “obtained index” 0 indicates that a total 
of 59 SPs of Ve having index 0 are found using the homotopy H^°\ and under the permutation symmetry they collectively 
represent two geometric configurations of the 6 atoms. No nondegenerate SPs of any other indices were found in these runs. 


N 

Vv at SP 

# neg. eig. 

# zero eig. 

Min. eig. 

5 

-3.0000000002858558 

0 

6 

10' 13 

5 

-4.0000000000004334 

0 

5 

io- 19 

5 

-6.0000000000157154 

0 

3 

10~ 14 

6 

-6.0000000007539667 

0 

6 

IQ- 13 

6 

-8.4806201568350339 

1 

3 

10 -12 

6 

-9.1038524157210112 

0 

3 

10- 14 

7 

-9.1038524162043597 

0 

6 

10- 13 

8 

-12.302927529932436 

0 

6 

10' 13 

8 

-14.811129645604781 

1 

3 

10- 13 

8 

-15.444733800924601 

1 

3 

10' 12 

8 

-15.533060054591045 

0 

3 

10- 15 

8 

-16.505384168030663 

0 

3 

10“ 15 

9 

-18.597348435556142 

1 

3 

10" 12 

9 

-18.778208174125020 

0 

3 

10- 14 

10 

-18.856826168132727 

0 

6 

10- 16 


Table II. A selection of numerically degenerate SPs of Vv 
obtained by the index-resolved fixed-point homotopy 0 for 
a range of N values. For each SP, “Vjv at the SP” is the 
value of Vv evaluated at the SP. “# neg. eig.” is the num¬ 
ber of negative eigenvalues of 7-L{Vn) at the SP. “# zero eig.” 
is the number of eigenvalues of 'H(Vn) at the SP that have 
magnitude less than 10~ 10 and do not correspond to the in¬ 
herent degree of freedom induced by rotation and translation. 
“Min. eig.” shows the approximate scale of the of the small¬ 
est magnitude of the eigenvalues counted in the column “ff 
zero eig.”. 


V. PRESERVATION OF INDEX UNDER RELAXED 
CONDITIONS 



Figure 2. The t value (horizontal) plotted against the arc- 
length (vertical) along a trajectory of the index-resolved fixed- 
point homotopy with N = 7 and m = 2. For each point (x, t), 
the index of x as an SP of V/ 2 ' is indicated by color: red/dark 
for index 2 and green/light for index 1. 

As discussed in 0 by the continuity of eigenvalues of 
H(Vf m) ) along the trajectory of x = a, the index must 
remain the same as long as no degenerate SP of V^ rn> is 
encountered (where 7^(V r / m ' 1 ) becomes singular). How¬ 
ever, our numerical experiments with the Lennard-Jones 
cluster (j jIV| ) suggest that the index can be preserved 
under much more relaxed conditions. In particular, we 
observed that index can he preserved even when turning 
points are encountered. Figure [2] shows the t -value plot 
together with the changes of indices along a curve defined 
by = 0 for the Lennard-Jones cluster potential 0) 
with N = 7 and target index m = 2. After four changes 
(from index 2 to 1 , then 2, 1 , and back to 2), the index 
comes back to the target index m = 2. 

Among our experiments summarized in Table |TJ such 
phenomenon where trajectories encounter some degener¬ 
ate SPs before reaching SPs of the target potential appear 
to be quite common. Indeed, we estimate that at least 
98% of the SPs counted in Table |T| for N = 10,..., 15 
(where all SPs obtained have indices agreeing with target 
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indices) were obtained by curves that encounter degen¬ 
erate SPs, suggesting that a much weaker condition may 
exist for the preservation of indices. 


VI. CONCLUSION 

Specializing for the potential energy landscape scenar¬ 
ios, we have proposed a novel homotopy continuation 
based method, called index-resolved fixed-point homo¬ 
topy, which can find SPs of specific index. It does not 
require a preexisting set of SPs as its starting points. 
The method can also find certain singular SPs of par¬ 
ticular index without much numerical difficulties unlike 
other Newton or quasi-Newton based methods. With our 
numerical experiments with the Lennard-Jones clusters, 
we have also demonstrated that the proposed homotopy 
continuation approach can target SPs of a specific in¬ 
dex with mild conditions on the potentials. This result 
may trigger further activities in such specialized homo¬ 
topy continuation based approaches among the relevant 
mathematics community. Future work will involves the 
rigorous analysis of this homotopy method as well as a 
systematic comparison of our proposed approach with 
various existing methods to find the SPs of a specific in¬ 
dex. 
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